%matplotlib nbagg
import numpy
numpy.random.seed(1)
import matplotlib.pyplot as plt
import matplotlib
import matplotlib.cm as cmx
import sdeint
import sys
sys.path.append('..')
from nlpl import *
from load_models import *
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import rc
plt.rc('text', usetex=True)
plt.rc('font', family='serif', size = 18)
p_max = 10
#p_max = 30
N = 10000
# model_name = 'nanopore'
# model_name = 'snanopore'
# model_name = 'lorenz'
model_name = 'slorenz'
# model_name = 'nanopore_iei'
# model_name = 'snanopore_iei'
# model_name = 'eeg'; roi_src = 10
# model_name = 'meg'; roi_src = 30
# model_name = 'gcp'; gcp_sensor_ind = 0 # Max is 26
nn_package = 'sklearn' # Deterministic, but slower.
# nn_package = 'pyflann' # Fast approximation, but non-deterministic.
# Determines the upper bound on the number of nearest neighbors to use
# for the computation of NLPL:
# pow_upperbound = 0.9
pow_upperbound = 0.5
# pow_upperbound = 0.75
# Determines the number of nearest neighbors to use in estimating
# the specific entropy rate and normalized q-step specific entropy
# rate.
pow_neighbors = 0.5
if model_name == 'eeg' or model_name == 'meg':
x, p_true, model_type = load_model_data(model_name, N, roi_src = roi_src)
is_multirealization = True
elif model_name == 'gcp':
x, p_true, model_type = load_model_data(model_name, N, gcp_sensor_ind = gcp_sensor_ind)
is_multirealization = False
else:
x, p_true, model_type = load_model_data(model_name, N)
is_multirealization = False
if model_name == 'nanopore':
x = x + (numpy.random.rand(x.shape[0]) - 0.5)*(1e1)
if 'iei' in model_name:
x = numpy.log(x)
p_opt, nlpl_opt, nlpl_by_p, er_knn, ler_knn = choose_model_order_nlpl(x, p_max, pow_upperbound = pow_upperbound, nn_package = nn_package, is_multirealization = is_multirealization, announce_stages = False, output_verbose = True)
print 'Chose p* = {}, giving ER(p*) = {}'.format(p_opt, er_knn)
plt.figure(figsize = (10, 5))
plt.plot(range(0, p_max + 1), nlpl_by_p)
plt.xlabel('Model Order $p$')
plt.ylabel('NLPL$(p)$')
plt.axvline(p_opt, color = 'r')
plt.tight_layout()
print 'Estimating specific entropy rate in-sample with q = 0...'
ser = estimate_ser_insample(x, ler_knn, p_opt = p_opt, pow_neighbors = pow_neighbors)
fig, ax = plt.subplots(2, 1, sharex = True, figsize = (10, 5))
ax[0].plot(x)
ax[0].set_ylabel('$X_{t}$')
ax[1].plot([numpy.nan]*p_opt + ser.tolist())
ax[1].set_ylabel('$\\widehat{{H}}^{{S}}(X_{{t-{}}}^{{t-1}})$'.format(p_opt))
ax[1].set_xlabel('Time Step $t$ (au)')
plt.tight_layout()
q = 5
print 'Estimating the predictive divergence in-sample with q = {}...'.format(q)
kl_q = estimate_normalized_qstep_insample(x, p_opt, q, pow_neighbors = pow_neighbors)
# Compare time series, SER, and normalized q-step SER for the full time series.
fig, ax = plt.subplots(3, 1, sharex = True, figsize = (10, 5))
ax[0].plot(x)
ax[0].set_ylabel('$X_{t}$')
ax[1].plot([numpy.nan]*p_opt + ser.tolist())
ax[1].set_ylabel('$\\widehat{{H}}^{{S}}(x_{{t-{}}}^{{t-1}})$'.format(p_opt))
ax[2].plot([numpy.nan]*p_opt + kl_q.tolist())
#ax[2].set_ylim([0, numpy.max(kl_q)])
ax[2].set_ylabel('$\\widehat{{\\widetilde{{H}}}}^{{S}}(X_{{t-{}}}^{{t-1}}; {})$'.format(p_opt, q))
ax[2].set_xlabel('Time Step $t$ (au)')
plt.tight_layout()
# Compare time series, SER, and normalized q-step SER for a snippet of the time series.
fig, ax = plt.subplots(3, 1, sharex = True, figsize = (10, 5))
ax[0].plot(x)
ax[0].set_ylabel('$X_{t}$')
ax[1].plot([numpy.nan]*p_opt + ser.tolist())
ax[1].set_ylabel('$\\widehat{{H}}^{{S}}(x_{{t-{}}}^{{t-1}})$'.format(p_opt))
ax[2].plot([numpy.nan]*p_opt + kl_q.tolist())
#ax[2].set_ylim([0, numpy.max(kl_q)])
ax[2].set_ylabel('$\\widehat{{\\widetilde{{H}}}}^{{S}}(X_{{t-{}}}^{{t-1}}; {})$'.format(p_opt, q))
ax[2].set_xlabel('Time Step $t$ (au)')
ax[2].set_xlim([0, 200])
plt.tight_layout()
ser_w_nan = [numpy.nan]*p_opt + ser.tolist()
kl_q_w_nan = [numpy.nan]*p_opt + kl_q.tolist() + [numpy.nan]*q
fig_range = [0, 200]
def uniformize(x):
a = numpy.nanmin(x); b = numpy.nanmax(x)
return (x - a)/float(b - a)
for measure_name, measure in zip(['Specific Entropy Rate', 'Normalized $q$-step Specific Entropy Rate'], [ser_w_nan, kl_q_w_nan]):
_uds = uniformize(x)
_udd = uniformize(measure)
_s = range(len(x))
_t = range(len(measure))
tspan = range(len(measure))
a = []
b = []
for s, t, u, v in zip(_s, _uds, _t, _udd):
a.append(s)
a.append(u)
a.append(None)
b.append(t)
b.append(v)
b.append(None)
plt.figure(figsize = (10, 5))
plt.plot(tspan, _uds, color = 'blue', linewidth = 2)
plt.plot(tspan, _udd, color = 'red', linewidth = 2)
plt.plot(a, b, color = 'green', linewidth = 0.25)
plt.xlim(fig_range)
plt.title(measure_name, fontsize = 12)
for measure_type in ['ser', 'kl_q']:
if measure_type == 'ser':
measure = ser_w_nan
cs = uniformize(measure)
elif measure_type == 'kl_q':
measure = kl_q_w_nan
cs = uniformize(measure)
fig = plt.figure()
colorsMap = 'jet'
cm = plt.get_cmap(colorsMap)
cNorm = matplotlib.colors.Normalize(vmin=numpy.nanmin(cs), vmax=numpy.nanmax(cs))
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)
ax = Axes3D(fig)
if measure_type == 'ser':
ax.scatter(x[:-3], x[1:-2], x[2:-1], s = 10, c=scalarMap.to_rgba(cs[3:]), lw = 0)
elif measure_type == 'kl_q':
ax.scatter(x[:-(3 + q)], x[1:-(2 + q)], x[2:-(1 + q)], s = 10, c=scalarMap.to_rgba(cs[3:-q]), lw = 0)
scalarMap.set_array(cs)
fig.colorbar(scalarMap)
# ax.set_xlim([-2, 2])
# ax.set_ylim([-2, 2])
# ax.set_zlim([-2, 2])
from matplotlib import animation, rc
from IPython.display import HTML
#X_anim = embed_ts(x, p_max = 2)[::2, :]
X_anim = embed_ts(x, p_max = 2)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot([X_anim[0, 0]], [X_anim[0, 1]], [X_anim[0, 2]])
# initialization function: plot the background of each frame
def init():
ax.clear()
ax.set_xlim([numpy.min(x), numpy.max(x)]); ax.set_ylim(numpy.min(x), numpy.max(x)); ax.set_zlim(numpy.min(x), numpy.max(x))
ax.plot([X_anim[0, 0]], [X_anim[0, 1]], [X_anim[0, 2]], 'b')
return (ax,)
def animate(i):
n = i + 3
ax.plot(X_anim[n:n+2, 0], X_anim[n:n+2, 1], X_anim[n:n+2, 2], 'b')
ax.scatter(X_anim[n, 0], X_anim[n, 1], X_anim[n, 2], color = 'b')
ax.scatter(X_anim[n+1, 0], X_anim[n+1, 1], X_anim[n+1, 2], color = 'r')
return (ax,)
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=200, interval=200, blit=False)
#HTML(anim.to_html5_video())
HTML(anim.to_jshtml())